clear;

opts.ar=[0.5 0.2 0.1]; 
opts.ma=0.6;
opts.mdlv=1;
dt=1;
t=0:dt:300;
mdl = arima('MA',opts.ma,'AR',opts.ar,'Variance',opts.mdlv,'Constant',0);
realizations=1;
x=simulate(mdl,length(t),'NumPaths',realizations);
[P s]=pmtmPH(x,dt,3,0);
P=log(P);

figure(1); clf; hold on;

pall=[1 1 2 2 3];
qall=[0 1 1 2 2];

for iter1=1:5,
    p=pall(iter1);
    q=qall(iter1);
    mdl2=arima(p,0,q);
    [mdlest mdlcov]=estimate(mdl2,x(:),'Display','off');
    %infer arma model
    mest=[mdlest.Constant mdlest.AR{1:end} mdlest.MA{1:end} mdlest.Variance]';
    
    %generate distribution of arma spectral fits
    num=1000;
    mestr=mvnrnd(mest,mdlcov,num);
    for iter=1:num,
        Pfit(iter,:)=log(arma_spec(s,dt,mestr(iter,2:1+p),mestr(iter,2+p:1+p+q),mestr(iter,end)));
    end;
    Pfito=log(arma_spec(s,dt,mestr(iter,2:1+p),mestr(iter,2+p:1+p+q),mestr(iter,end)));
    
    plot(s,P,'k','linewidth',1); 
    plot(s,Pfito,'b'); 
    plot(s,prctile(Pfit,50),'r'); 
    if iter1==3,
        plot(s,prctile(Pfit,95),'r--'); 
        plot(s,prctile(Pfit,05),'r--');
    end;
    axis tight;
    drawnow;
end;


return;
%combine two distribution through convolution

n=10^6;
x1=randn(n,1);
x2=randn(n,1).^2+randn(n,1).^2;
x=x1+x2;

yi=-25:0.001:25;
y1=normpdf(yi,1,1);
y2=chi2pdf(yi,2);
yc=conv(y2,y1,'same'); 
yc=yc/(sum(yc)*0.001);

figure(2); clf; hold on;
ksdensity(x,yi);
%plot(yi,y1,'r');
%plot(yi,y2,'g'); 
plot(yi-1,yc,'m');

